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SUMMARY 

A continued fraction representation for Theodorsen's circulation function 
is derived. This continued fraction converges to the circulation function 
everywhere except on the branch cut. It can be used to compute the function 
except when the argument is small. When converted to pole-residue form the 
continued fraction greatly facilitates the evaluation of integrals containing 
the circulation function. 


INTRODUCTION 

Theordorsen 1 s circulation function (ref. 1) relates lift to downwash 
in unsteady incompressible potential flow. The function can be expressed as 
the ratio of two continguous confluent hypergeometric functions and hence has 
a continued fraction representation derivable from the continued fraction of 
Gauss. 

This continued fraction can be truncated to give rational approximations 
to the circulation function. These approximations are useful in control 
theory because their poles and zeroes are easily computed. These approxima- 
tions can also be inverse Laplace transformed to give accurate approximations 
to Wagner's function. 


SYMBOLS 

A recursion coefficient matrix 

A 

n 

A $,k' B Z 


a k' b k' c k 
B 

n 

b 

n 

C(-is) 


truncated continued fraction approximation to C(-is) 
least square approximation to C(03) 


C n ( - is) 
c (00) 


truncated continued fraction numerator 
coefficients of least squares simultaneous equations 
continued fraction coefficients (numerator) 
polynomial recursion coefficients 
truncated continued fraction denominator 
continued fraction coefficients (denominator) 

Theodor sen's circulation function 


//&- 32 . 3*30 ^ 


E 2„ (t) 


_F (a, b; ; z) 
z o 

F (-is) 

G(-is) 


V Y v 


n 


n' 


n 


n 


n 


R 


\ <x » 


s 

s k 

I 

s k 

t 

x 

x 

At 


u, t v, 
k k 


diagonal elements of A-matrix 

correction integral used when evaluating Wagner’s function 
subdiagonal and superdiagonal A-matrix elements 
confluent hypergeometric function 

even part of C(-is) (real part of C(-is) if -is is real) 
[c(-is) - F(-is)]/(i2) 

modified Bessel function of the first kind 

/T 

Bessel function of first and second kind 
truncation order 

number of residues that contribute significantly to C* (-is) 

2n 

even truncation numerator 
odd truncation numerator 
even truncation denominator 
odd truncation denominator 

ratio of two contiguous confluent hypergeometric functions 

any polynomial of degree k satisfying a three term 
recursion relation 

complex argument of circulation function = 0 + ito 

poles of (-is) 

2n 

zeros of numerator in (-is) 

2n 

time 

polynomial column vector 
4s 

result spacing in FFT quadrature 

coefficients used in least squares approximation 
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4> (t) Wagner's function 

to a value of to such that |c„ (to) - C(to) I < e if to > to 

o zn o 

Aw step size used in FFT quadrature 

asymptotic to 

* approximately equal 

X eigenvalue of A (= -4s) 

V arbitrary order of Bessel function 

0 real part of s 

GO imaginary part of s 

THE CONTINUED FRACTION 

Theodorsen's circulation function can be expressed 


C (-is) = K. (s) / (K (s) + K. (s) ) 
1 o 1 


(1) 


The region of aerodynamic interest in the complex s-plane lies on or near the 
positive imaginary axis. The Bessel function K^(s) is expressible as a 
confluent hypergeometric function (eqs. 13.6.21 and 13.1.10.2 of ref. 2). 



(2) 


Replacing the Bessel functions in C = 1 - K q /(K o + K ) by confluent hyper- 
geometric functions and using Gauss' relation for contiguous functions 
(eq. 15.2.14 or 13.4.17-18 of ref. 2) to combine the numerator functions gives 


C (-is) = 1 - h 




(3) 


The ratio of two contiguous confluent hypergoemetric functions 


R = 0 F (a, b; ; z) / F (a, b - 1?; z) (4) 

2 o 2 o 

has a very simple continued fraction representation (the confluent form of the 
continued fraction of Gauss, see chapter XVIII of ref. 3). It is 
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R = 


1 

az 


1 - 


1 - 


bz 


1 - 


(a + l)z 


1 - 


(b + l)z 


1 - 


(a + 2)z 

1 - . . 


(5) 


so 


C (-is) = 


b + 
o 


b^ + 


b 2 + b 7T 


where b Q =1, a 1 = -h, b 1 = 1 


(6) 

(7) 


and 


a 2n = 2n " 


2n 


4s, 


n 


1 2 

r ^ t • • • 


(8) 


a„ 

2n+l 


2n 


1, 


b 2n+l 


1 


That is, 


C (-is) 


1 - 


h 


1 + 


4s + 


1 + 


4s + 


(9) 


This continued fraction converges to C(-is) in the entire complex s-plane 
cut on the negative real axis. 
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RATIONAL APPROXIMATIONS 


Let C n (-is) represent the continued fraction of equation (6) and (7) 
truncated by discarding all terms beyond a n /b n . That is, by setting a n +i 
to zero. This can be expressed as a rational function 


C = A /B (10) 

n n n 


where A n , B n are polynomials in 4s computed by the usual forward recur- 
sion formula (ref. 3) for continued fractions, namely 


A i = 1; B = 0; A = b ; B = 1 
-1 -1 o o o 


\ ~ b k\-l + a k\-2 ; B k b k B k-l + a k B k-2' k 1,2 "“ 


(ID 


Since the even and odd terms of the fraction have different forms (eq. (8) ) , 
it is convenient to separate the even and odd subscripts in recursion (11) . 
Let 


C 0 (-is) = * 5 P (4s) /Q (4s) 
zn n n 


C, (-is) = hP (4s)/Q (4s) 
2n+l n n 


( 12 ) 

(13) 


where 


and are 


2S 2„' 

2n = 

B 2n 



2A 2„ + 1 ! 

2n 

= B 2n+1 


(14) 

formulas 

for 

p n» Qn' p n 

and Q n 

are derived from equation (11) 

:) = 2; 


Q 0 ( x) 

= 1 


:) = x + 

2; 

Q, M 

= X + 1 

(15) 
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P k+l (x) = (x + 4k > p k < x ) - < 2k - 1)2p k _i (x ) 

Q k+1 (x) = (x + 4k)Q k (x) - (2k - l^Q^x) (16) 

P (x) = 1 Q (x) = 1 

o o 

P-^x) = x + 3 Q x (x) = x + 2 (17) 

P. .. (x) = (x + 4k + 2) IP (x) - (4k 2 - 1)P. (x) 

k+1 k k-1 

Q k+ i ( x) = (x + 4k + 2)Q k (x) - (4k 2 - 1)F (x) (18) 

It can be seen that P n , Q n , P n , and Q n are all polynomials in x = 4s 

of degree n for n > 0. The first few polynomials of each set are 

Q=1 Q = x + 1 

*o *1 

Q_ = x 2 + 5x + 3? Q- = x^ + 13x^ + 34x + 15 

2 3 

Q. = x 4 + 25x 3 + 165x 2 + 298 x + 105 

4 

= x 5 + 41x 4 + 516x 3 + 2301x 2 + 3207x + 945 (19) 

5 

P = 2; P = x + 2 

o 1 

2 3 2 

P^ = x + 6x + 6; P^ = x + 14x + 45x + 30 

P = X 4 + 26x 3 + 188x 2 + 420x + 210 

4 

P 5 = X 5 + 42x 4 + 555x 3 + 2742x 4 + 4725x + 1890 (20) 
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Q 0 - 1; Q 1 = x + 2 

Q = X 2 + 8x + 9; Q. = X 3 + 18x 2 + 74x + 60 

2 d 

Q = x 4 + 32x 3 + 291x 2 + 216x + 525 
*4 

O = x 3 + 50x 4 + 804x 3 + 4920x 2 + 10551x + 5670 (21) 

5 


P = 1; P = x + 3 

O -L 

p = x 2 + 9x + 15; P = X 3 + 19x 2 + 90x + 105 

2 . j 

P = X 4 + 33x 3 + 321X 3 + 1050x + 945 
4 

P = x 5 + 51x 4 + 852X 3 + 5631X 2 + 14175x + 10395 (22) 

Inspection of the first new polynomials and the recursion formulas shows that 

Q (0) = (2n - 1) !! 

P (0) = 2(2n - 1) ! ! 
n 

Q (0) = (n + 1) (2n - 1) ! ! 

P (0) = (2n + 1) ! ! < 23 ) 

n 


c 2n (°) = 1 


and 


C 2n+1 (0> 1 2n + 
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The even numbered convergents C2 n give the correct value of C(0) while 
the odd convergents C2 n +1 merely approach the correct value. Because of 
this, and because they have a slightly simpler eigenvalue matrix, the even 
convergents are much more convenient to use. 

The even convergents are the diagonal elements, and the odd convergents 
are subdiagonal elements of a Pade matrix defined by setting its first column 
to the convergents of the asymptotic series for C(-is) and setting its first 
row to those of the asymptotic series for 1/C (-is). Reference 4 tabulates 
the first few diagonal Pade elements. However, the expression for CQ(-is) 
in reference 4 is incorrect. Fortunately reference 4 makes no further use 
of Cg(-is). 


POLES AND ZEROS 

All the poles and zeros of C2n(“ is ) li e on the negative real axis. 
If they are numbered in order of distance from the origin they satisfy the 
inequality 


_oo < s‘ < s < s ' ' . . . < s ’ < s 0 < s 1 < S -1 < 0 

n n n-1 2 2 1 -L 

where s^ are the poles (zeros of Q n ) and s^- are the zeros (zeros of P n ) . 

It is not practical to compute the poles of C 2 n (-is) by solving 
Q n (x) =0 as a polynomial equation if n is large because Q n (x) overflows 
the computer if x is barely outside of the convex set containing all of the 
roots. Instead the recursion formula (16) is used to construct matrices 
whose eigenvalues are the poles and zeros of C2 n (-is). 

Suppose the polynomials Rj^(x) are each of degree k for k 0 and 
satisfy the three term recursion relation 


"kVi - (x + b k )E k + 'A-i - 0 1241 

If R^i / 0 then a Q and b Q should be redefined so R_i does equal 0. If 
R n (x) = 0, then equation (24) can be written in matrix form as 


8 


— — 


— — 


— — 

b -a 0 0 0 0 


R 


R 

o o 


o 


o 

-c, b, -a, 00 0 


R, 


R, 

111 


1 


1 

0 

o 

1 

03 

1 

( 

A 

< 

o 

o 




R_ 

2 2 3 


2 


2 

o 

A 

< 

0 

1 

o 

o 


R_ 

= -X 

R_ 

J 3 


3 


3 

0 





. 





0 0 0 0 0 -c _ b _ 


R 


R «, 

n-1 n-1 


n-1 


n-1 

— 








( 25 ) 


This is a matrix eigenvalue problem 


AX = Ax 


( 26 ) 


where X = -x 


( 27 ) 


X = col (R , R, , . . . R . ) 
o 1 n-1 


( 28 ) 


and A is a tridiagonal matrix whose kth row is 


(O' *•' * c k-i' b k-i' ~ a k-i' •*- 0) 

The matrix A is made symmetric by replacing R^ by where 


( 29 ) 


Y k /Y k-1 >/c k /a k-l 


( 30 ) 
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Then AX = Ax where now 


( 31 ) 


X = col (R , R , ... R . ) 

o 1 n-1 

and the kth row of A is 


( 32 ) 


(O' **•' V V e k+i' ••• 0) 


( 33 ) 


where 


d k b k-l 


( 34 ) 


as before and 


e k - /c k-l\-2 


( 35 ) 


To compute the poles of C^^C-is) let = Q . Then 


A = 


1 

-1 

0 

0 


-1 

4 

-3 

0 


0 

-3 

8 

-5 


0 

0 

-5 

12 


0 

0 

0 

-7 


( 36 ) 


and s k = ~ hX k‘ 


( 37 ) 
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That is 


for 


d i = 1 


e 2 " 1 


d 2 = 4 


k = 3 to n 


d k = d k -i + 4 


e. = e - 2 
k k-1 


next k 


(38) 


Similarly the zeros s k of C 2n (-is) are obtained by letting = P k « 
Then 


A = 


2 

-/2 

0 

0 


-/2 

4 

-3 

0 


0 

-3 

8 

-5 


0 

0 

-5 

12 


0 

0 

0 

-7 


(39) 


s k = " u k 


(40) 


11 



That is 



for k = 4 to n 

*k " d k-l + 4 


next k 


-\ 


> 


J 


(41) 


Both A matrices are tridiagonal f symmetric, and positive definite. Eigen- 
values can easily be computed even for very large n because the nonzero matrix 
elements are easy to compute and do not vary widely in magnitude. Table I was 
computed using procedure tg£l on page 232 of reference 5. It lists the 

poles s and zeros s' of C (-is) for n = 1, 2, 4, and 8. 
k k zn 

POLES AND RESIDUES 

The expression for C 2 n (-is) can be written 


c 2 „ ( -is> 


P (4s) - Q (4s) 

JL n ■ n 

2 2Q (4s) 
n 


(42) 


P n " Q n is of lower degree than Q n and the zeros of Q n are all distinct 
so, using partial fractions, one obtains 


C 2n (-is) 



(43) 


The coefficients of the partial fraction r^ are the residues of the poles 
s^. They can be computed from either 
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V 4 s k> 


(44) 


or 




r, = 


2 <s k - <> n 


S,_ - S r 


k £4 s k 

JWk 


- s„ 


(45) 


The residues approach zero rapidly as k increases. This is because of the 
factor Sk - in equation (45) . Except for the upper left hand corner 
elements d-^ and e 2 the eigenvalue matrices for and s £ are identical. 

These corner elements have little effect on the higher eigenvalues because the 
matrices are diagonally dominated. The fact that r^ -*• 0 rapidly as k 
increases is important because it means that when n is large the sum (43) 
can be truncated as n' where n' < n and 


IV s k! < e 


(46) 


for all k > n'. Then 


C 2n ( - iS) 




(47) 


The number of terms in the sum (47) increases much more slowly than n. 
Table II , which is a continuation of table I, lists the poles, zeros, and 
residues for n = 16, 32, 64 and 128, It also lists n 1 based on 
e = io~10. Note how slowly n* increases with n. 


APPLICATIONS 

Applications of the continued fraction representation of C(~is) include 
its use to evaluate the function and its use to represent or to evaluate 
integrals (particularly infinite limit integrals) containing C(-is) . 

The error contours of C2n(~i s ) resemble a family of parabolas containing 
the negative real axis and with a common focus at the origin. That is, if 
s = a + ico (a,to real) , then the error contours are approximated by the family 
of parabolas 
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(48) 


2 

2cj a + to = to ' 

o o 

C0 o , the intercept of the parabola with the co-axis, is a function of e, the 
error tolerance , and of n. Given n and £ if to Q has been computed then 
|c 2 n (~is) - C(-is) | < £ for all s = a + ito for which 


2 2 

(0 + 2(0 a > to (49) 

o o 

It is hard to compute (0 Q from £ and n so the usual procedure is to 
choose (0 o and then compute 

e = [c_ (oj ) - c(u> ) | (50) 

r 2n o o 1 

-12 

For example, if to Q = 2 and n = 8, then £ = 0.3 x 10 and the error in 
C 2 n is less than £ for all real (0 > 2 and for all complex s = 0 + ito 
for which 4a + co 2 > 4. 

The continued fraction is the most efficient way available to compute 
C (co) if (0 _> 2 or if complex s = 0 + i(0 satisfies 40 + ( 0 2 4. It should 

be used in pole-residue form truncated at n* for small co. If co > 20 
(or if 400 + co 2 > 400) it should be used as a truncated asymptotic series. 

The asymptotic series for C(-is) is 


C (-is) 



2 7 

2 + 3 

(4sr (4sr 


38 

(4s) 4 


286 2756 32299 

- + ~ 

(4s) (4s) (4s) 


(51) 


as s 00 for |arg(s) | < — . 

The asymptotic series is obtained from equation (9) by repeated division. 

The series (51) diverges for all s. The asymptotic series, unlike the 
continued fraction, can be used to approximate C(-is) on the branch cut, 
arg(s) = ±7T, if |s| is sufficiently large. 

The other use of the continued fraction mentioned at the beginning of this 
section is to facilitate evaluation of integrals containing C(-is). For this 


14 , 



application the pole-residue form, equation (47) , is used. Two examples are 
presented. One is the evaluation of Wagner's function <j)(t). The other is 
the evaluation of some integrals that occur when approximating C(-is) using 
least squares. 


Wagner's Function 

Wagner's function (f>(t), is the inverse Laplace transform of C(-is)/s. 


4> (t) 


1 

i2TT 


/•c+i«> 

% 'c-ioo 


st C (-is) 

e 

s 


ds 


(52) 


Each rational approximation to C(-is) 


C 2n ( ' is) 


n*<n 



(53) 


has an associated exponential approximation to cf) (t) 


n'<n 

r. 


*2n lt> " 1 + E 

k=l 


k e s k t 


(54) 


obtained by substituting equation (53) into equation (52) and performing the 
indicated integral transformation. The 1 appearing in equation (54) is 
computed from 


C 


2n 


( 0 ) 


n 


1 

2 





(55) 


This is only true for the diagonal Pade elements C 2 n . For the subdiagonal 
Pade elements 
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C 2n+1 


( 0 ) 


1 

2 





1 

2n + 2 


(56) 


n'<n 


W 4 ’ - 


1 - 


2n + 2 


k=l 



(57) 


The odd approximants to (j)(t) do not give the correct limit at t=°°. The 
approximation (54) can be used to compute cj> (t) accurately, even for small n, 
if t is sufficiently small. For large t equation (54) has too strong an 
exponential decay. 

Equation (54) illustrates the use of the pole-residue form of the 
continued fraction for C(-is) to replace a numerical integration by a closed 
form integration. It can also be used to simplify a numerical integration. 

To integrate equation (52) the path of integration (c-i°°, c+i<*>) must be 
fixed. The only restriction is that the path be to the right of the branch 
point at s=0 and have a nonpositive real part at the two ends. Two paths 
are very convenient for numerical integration. One is the imaginary axis as 
shown in figure 1. The other is the branch cut as shown in figure 2. 

If path 1 is used and symmetry of the integrand is considered (see 
sections 5 through *7 of ref. 6 for details) one obtains the following integral 
representation for cj)(t) 


>.00 


<Mt) =- 


sin cot 
0) 


F (CO) dlO 


(58) 


If path 2 is used and symmetry of the integrand is considered one obtains 


4> (t) = l - 


“2 -st , 
s e ds 


I (K o - K 1> 2 + ,2|I „ + V 2 


(59) 


The continued fraction for C(-is) cannot be used to help evaluate 
equation (59) because the continued fraction diverges along the branch cut. 
However, the continued fraction can provide . considerable help in evaluating 
equation (58) as wil'l be shown. 
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The function F (to) appearing in equation (58) is usually defined to be 
the real part of C(to). However, if complex arguments are permitted, it is 
more convenient to define it to be the even part of C(-is). 


F (-is) = ^{c (-is) + C (is) ] 


(60) 


The function F(-is) usually occurs in association with the function 
G(-is) defined 


G (-is) 


utc(-is) 


C (is) ] 


(61) 


Equation (53) furnishes rational approximations to F and G 


F_ (to) 
2n 


n'<n 

S 1 + 
k=l k 

n'<n 


(62) 


03 


G (03) = 03 
2n 


£ 


(63) 


k=l S k 


03 


If the F(03) in equation (58) is expressed 


F (to) = F (tO) + [F(tO) - F ((0)] 
2n zn 


(64) 


and the first F^ n is integrated in closed form one obtains 


4>(t) =<f> 2n (t) + E 2n (t) 


(65) 


where (J> 2 n (t) is given by equation (54) and the error or correction term 
E 2n (tl is 


E 2„ lt> - - 


»r 


sin tot 

to 


[f (to) - F 2n (u)] dto 


( 66 ) 
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Equation (66) is much easier to integrate numerically than equation (58) . 
This is because it has finite limits and because a large relative error can 
be tolerated. Infinite limit oscillatory integrals are notoriously difficult 
to compute. Equation (66) has a finite upper limit because F - F 2 n is 
essentially zero for u) _> A large relative error can be tolerated because 

|F - F 2n | « F even for w < w 0 . 

If equation (66) is to be integrated for a single value of t a 
sophisticated integration technique such as Legendre-Gauss or Romberg quadrature 
can be used. However, if a large number of values of t are used, then 
equation (66) should be evaluated as a trapezoidal sum using a fast Fourier 
transform. For a quadrature order m let 



0 ) = .kAio, 


k = 0 to m 


(67) 


t = t + £At, 
o 


£ = 0 to m 


The FFT formalism requires that 


AojAt = 27T/m 


( 68 ) 


so 


a 2TT 

At = — 

03 

o 


(69) 


The trapezoidal sum for E^ n is 


E 2„(to + «it) 


m 

. M «« 

— Aw / sin kAw (t + &At) 

TT / J O 

k =0 


F(kAw) - F„ (kAw) 
zn 

kAw 


(70) 
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or 


E ?n (t + ^ SL) ~ -£m 
2n o U) 

o 


m-l 

2 X ' -i2ir£k/m , 
?2^ e f k 

k=0 


for Z - 0 to m-l 


where 


f = 0 
o 


and 


■ W O t O v 

-i k 

m 


f k = 




for k = 1 to m-l. 


(71) 


(72) 


By using a fast Fourier transform to evaluate the sum in equation (71) , 
it is possible to evaluate E2n( t Q * 27T / (j0 o ^ for m va l ues of & using 
only log 2 (m) times as much computing effort as would be required for one 
value of Z. However , only the E 2 n (t Q + 2 tt/cj0 o Z) for Z < m/4 are reason- 
able approximations to the integral (66) . The term t Q in the argument of 
E 2 n is to permit interpolating between values of Mt. It should not be 
larger than 2 tt/oj 0 or aliasing can occur. Thus, the largest value of t for 
which 4) (t) can be computed is 


t 


7T m 
2w 


o 


(73) 


This can be increased either by increasing the quadrature order m or by 
increasing the exponential approximation order n, thereby decreasing 0) q . 

The circulation function F(w) in equation (72) can be computed using 
Bessel functions 


F(W) = 


J 1 (J 1 * V * V Y 1 - J o> 

( J 1 + v 2 + (X 1 - v 2 


(74) 


A very convenient way to compute or Y^ for a large number of equispaced 

arguments 
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k = 1 to m - 1 


(75) 
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is to evaluate Jv and Yv accurately at the two ends, to = co 0 /m and 
to = (0 o - (0 o /m and then to approximate J V ' Y V at all the intermediate 
points (k = 2 to m - 2) by solving the Bessel equation as a finite difference 
boundary value problem. This permits computing J v , Y v at the intermediate 
points with less computing effort than is required for an elementary function 
such as a square root. It has the additional advantage that the accuracy 
with which the intermediate Bessel functions are computed increases as the 
quadrature order m is increased. 


Equations (65) , (53) , and (71) can be used to compute <}>(t) for any 

reasonable value of t. For unreasonable values of t such as 10^ or to 
estimate the behavior of 4> (t) as t 00 the other integral representation, 
equation (59), 'should be used. It can be evaluated as a Laguerre-Gauss 
quadrature. When setting up a Laguerre-Gauss quadrature the exponential 
weight function should represent the behavior of the integrand at infinity 
rather than merely being an. easily identifiable exponential factor. That is, 

- ( 2+t) s -st 

in equation (59) the exponential weight function is e , not e 

-2s 

The factor e comes from the asymptotic behavior of I q and 1^. Thus 


4>(t) 


1 



- (2+t) s 
e 


f (s) 


ds 


(76) 


where 


f (s) 


2s 


S 2 [(K - K ) 2 + IT 2 (I + I ) 2 ] 

o 1 o 1 


(77) 


Letting s = x/(2 + t) gives 


(t) = l 



( 78 ) 


which can be integrated numerically to give 
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( 79 ) 


m 

* (t > = 1 - ‘(itt) 

k=l 


where x^, w k are the Laguerre-Gauss abscissas and weights. As t -t 00 
the integral in (78) approaches 



e X f(0) dx 


1 


(80) 


so a very Crude approximation to cf) (t) is 


<J>(t) 


1 - 


1 

2 4* t 


(81) 


A slightly better approximation is obtained by retaining the next term in 

the expansion of f(s) about the origin and integrating in closed form. This 

gives 


4>(t) 


1 


1 

2 + t 



2 

2 4* t 


a n (4 4- 2t) 


(82) 


Equations (81) or (82) show that <f>(t) approaches 1 like 1/t as t 
increases rather than exponentially as indicated in equation (54) . 


A Least Squares Approximation to C(-is) 

This section describes the use of the pole-residue form of the continued 
fraction for C(-is) to evaluate some integrals that occur when generating 
a least squares rational approximation to C(-is). 

The continued fraction for C(-is) is an expansion about s = 00 4- iO 
and hence is very accurate for large 0 and |to| (where s = a 4* ia)) and 
very uneconomical when O and co are both small. It has been shown 
(refs. 4 and 7) that economical rational approximations valid for 0 <_ 00 00 

and 0-0 can be generated by least squares. 

There are several ways of generating a least squares rational approximation 
to C(-is). The method described here has a moderately complicated derivation 
but is computationally very simple. 
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Let 


C(W) 


$ + Y 

2 v + io) 

k+1 


(83) 


where the constants and are chosen so C(o)) approximates C(o)) 

over the entire positive real co-axis. One way of computing the constants u k 
and v^ is to minimize the error 


0 ) I C (03) - C( 0 ))| 2 do) (84) 


The expression for C(co), equation (83) , is linear in the u^ and the nonlinear 
in the v^. Hence, it is very easy to minimize E with respect to the u^ 
and very hard to minimize E with respect to the v^. For this reason the 
Vfc are preassigned and the minimization is performed with respect to the u^ 
only. The v ^ are chosen to be -s^ for n* = m if no better choice is 
available. The error expression E, equation (84) , contains a weight function 

_i^ 

0) 2 . This is included to origin weight the error. It is needed because C (co) 
has a logarithmic branch point at co = 0 and hence is hard to approximate 
for small co. All that has to be done to compute the u^ is to set 9 e/3u^ 
to zero and solve the resulting system of linear simultaneous equations. This 
is performed as follows. 

Let 



C ( 03 ) = F (to) + iG(cO) 


(85) 


where 


F (to) 


m 


1 

2 


V=1 


u k V k 



( 86 ) 


G (to) 


m 


-“E 

k=l 


+ CO 


(87) 
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Then 


E 



of 1 * [ (F - F) 2 + (G - G) 2 ] 


dco 


(88) 


and 


1 _3e_ 

2 3„ t 




(G - G) 



dco 


0 


for Jl = 1 to n 


(89) 


This gives the system of equation 


m 


E 

k=l 


\k U k = (il = 1 to m) 


(90) 


where 


A £k “ 


tt/2 J 


0) 


v £ v k + w 


, 2 2 . . 2 2 . 

(v^ + co ) (v + co ) 


da) 


(91) 






0) 


V £ (F(C0) - -) - COG (CO) 

2 2 
V^ + CO 


dco 


(92) 


The factor 1/ (tt/ 2) was appended to both A^ and because it makes 

some subsequent arithmetic easier. 

The integral A^ can be evaluated in closed form. It is 


ilk 


v il + V k 


1 - + ^- 


(93) 
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The integral B£ has to be computed numerically. A very easy way to perform 
this numerical integration is to replace C = F + iG by C 2 n , the pole-residue 
form of a truncated continued fraction, with n chosen large enough so that 
the error C2n (w) - C (to) is negligible compared to the error C (to) - C (to) . 

Then 



(94) 


After A^j. and have been computed, equation (90) can be solved for u k . 

The above procedure minimized E with respect to u k for preassigned v k . 
If one wants to minimize with respect to v k , then substitute equation (90) 
into (88) to get the penalty function 


E = E 

o 



k=l 


(95) 


and perform a nonlinear optimization using some technique such as the Davidon- 
Fletcher-Powe 1 1 algorithm (ref. 8). The constant E Q in equation (95) is 
given by 


/•°° 



It can be computed using C 2n (-is) the same way B^ was computed. However, 
since it is a constant, it is not needed for the optimization and can be 
set to zero in equation (95) . 


CONCLUDING REMARKS 

The principal result of the investigation is the fact that Theodorsen's 
circulation function has a continued fraction representation with a particularly 
simple coefficient pattern, namely the consecutive odd integers. Although 
this continued fraction converges extremely slowly it still furnishes an 
economical way to compute the circulation function. The reason for this is 
that when converted to pole-residue form the terms containing distant poles 
can be discarded. For example, retaining only the 120 poles nearest to the 
origin from a 2048 term pole-residue form of the continued fraction introduces 
an additional error of only 10~10 to the approximation to C(-is). 
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This investigation also furnished some information about the singularities 
of the circulation function. The Bessel functions that comprise C(-is) , 
namely K^s) and K±(s) , have no singularities other than logarithmic 
branch points at s = 0 and s = °°. The only singularities that C(-is) can 
possess are these two branch points and, possibly, poles at the zeros of the 
denominator in C(-is), namely K^s) + Kj^s). Since C2 n (“is) converges to 
C(-is) everywhere except when arg(s) = ±TT , these zeros, if they exist at 
all, must all lie on the negative real axis. On the negative real axis the 
Bessel functions K Q and can be expressed as functions of positive 

real argument by analytic continuation (see eq. 9.6.31 of ref. 2). The 
real part of this analytic continuation is the same for all sheets of the 
Riemann surface and is Kq(x) - (x) where x = -s is real and positive. 

Inspection of figure 9.8 of reference 2, and of the asymptotic expression for 
K Q and K i, shows that K Q (x) - K-^(x) has no positive real zeros. Hence 
C(-is) has no singularities whatsoever except for logarithmic branch points 
at s = 0 and s = 00 . 

The continued fraction representation of C(-is) furnishes a very 
convenient way to compute Wagner's function from its definition as an inverse 
Laplace transformation. This leads either to an exponential approximation 
to Wagner's function or to an exponential approximation with a numerically 
integrated correction term. The latter can be used to compute <f> (t) (but 
not 1 - (j)(t)) to full register accuracy. 

Another application of the continued fraction that was discussed 
involved its use to evaluate some integrals that occur when approximating 
C(-is) by least squares. Low order approximations obtained using least 
squares are much more accurate over the frequency range of aerodynamic interest 
than the same order truncations of the continued fraction. However, high 
order approximations (those with over twenty terms) should be obtained by 
truncating the continued fraction because of numerical problems associated 
with the least squares process. 


25 



REFERENCES 


1* Theodor sen, T. : General Theory of Aerodynamic Instability and the 

Mechanism of Flutter, NACA Report 496, 1935. 

2. Abramowitz, M. , and Stegun, I. A., eds: Handbook of Mathematical Functions 

with Formulas, Graphs, and Mathematical Tables. NBS Appl. Math Ser. 55, 
U.S . Dep. Com., June 1964. 

3. Wall, H. S.: Analytical Theory of Continued Fractions, D. Van Nostrand, 

Inc., 1948. 

4. Vepa, R. : On the Use of Pade Approximants to Represent Unsteady Aerodynamic 

Loads for Arbitrarily Small Motions of Wings, AIAA Paper 76-17, 

AIAA 14th Aerospace Sciences Meeting, Washington, DC, Jan. 1976. 

5. Wilkinson, J. H. , and Reinsch, C. : Linear Algebra. Springer - Verlag, 

New York, -1971. 

6. Bisplinghof, R. L. , Ashley, H. , and Halfman, R. L. : Aeroelasticity, 

Addison-Wesley Inc., 1955. 

7. Dowell, E. H. : Unsteady Aerodynamics in the Time Domain, NASA TM-81844, 1980. 

8. Davidon, W. C. : Variance Algorithm for Minimization, Computer Journal, 

Vol. 10, No. 4, February 1968, pp. 406-411. 


26 




Figure 1 
Imaginary Axis 



Figure 2 
Branch Cut 
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TABLE I 


k 

-s k poles 

1 

zeros 

residues 

n = 1 


0.2500000000 

0.5000000000 

0.1250000000 ! 

n = 2 

i 

0.1743060906 

0.3169872981 

0.0798343811 

2 

1.0756039094 

1.1830127019 

0.0451656189 


n = 4 


1 

0.1156902163 

0.1742812056 j 

0.0414478723 

2 

0.5376897142 

0.7031255226 

0.0728261277 

3 

1.6861213269 

1.7116460071 

0.0105184411 

! 4 

i 

! 

3.9104987426 

3.9109472647 

0.0002075590 


n = 8 


i .1 

• 0.0708292313 

0.0895316590 

0.0162792778 

; 2 

j 0.3014411613 

0.4253867638 

0.0677288394 

3 

} 0.8103741892 

0.8990325533 

0.0333426883 

4 

; 1.6931876604 

1.7102975420 

0.0069379972 

5 

! 2.9738355808 

2.9753615342 

0.0006840380 

i 6 

4.7259501453 

4.7260072766 

0.0000268518 

7 

! 7.1126863436 

7.1126869821 

0.0000003071 

8 

10.5616956880 

10.5616956889 

0.0000000005 
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TABLE II. 


k 

-s^ poles 

“S. 1 zeros 
k 

r. residues 
k 


n = 16 

n' = 12 


i 

.0399290405 

.0450767705 

.0049890236 

2 

.1770117985 

.2309790353 

.0410711992 

3 

.4213665832 

.5226192420 

.0460005529 

4 

.8319542602 

.8955477962 

.0226896684 

5 

1.4174350819 

1.4379146184 

.0078656754 

6 

2.1767569894 

2.1814331978 

.0019842800 

7 

3.1151752001 

3.1159585681 

.0003530589. 

a 

4.2433929936 

4.2434850891 

.0000429141 

9 

5.5769738664 

5.5769811202 

.0000034466 

10 

7.1372634130 

7.1372637778 

.0000001755 

11 

8.9538356848 

8.9538356959 

.0000000054 

12 

11.0691194848 

11.0691194850 

.0000000001 


n = 32 

3 

It 

H 


1 

.0212539667 

.0225775151 

.0013308165 

2 

.1014880467 

.1181364893 

.0161066554 

3 

.2313703318 

.2837315435 

.0367590439 

4 

.4281625584 

.5028891998 

.0331916084 

5 

.7083237338 

.7663017562 

.0202106302 

6 

1.0728150641 

1.1020615789 

.0103265831 

7 

1.5200677586 

• 1.5318978031 

.0045717209 

8 

2.0496021188 

2.0537719222 

.0017407160 

9 

2.6619849235 

2.6632677368 

.0005642149 

10 

3.3585132712 

3.3588528353 

.0001544540 

11 

4.1409982347 

4.1410745802 

.0000354967 

12 

5.0116623678 

5.0116768049 

.0000068133 

13 

5.9731199112 

5.9731221892 

.0000010866 

14 

7.0284055993 

7.0284058970 

.0000001431 

15 

8.1810313620 

8.1810313940 

.0000000155 

16 

9.4350639127 

9.4350639155 

.0000000014 

17 

10.7952245954 

10.7952245956 

.0000000001 


n = 64 

n' = 23 


1 

.0109615336 

.0112936532 

.0003363445 

2 

.0550659069 

.0594029977 

.0045647906 

3 

.1282525325 

.1451687092 

.0165007629 
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TABLE II.- Continued 



n = 64 

n' = 23 


4 

.2294179265 

.2662436189 

.0271832071 

5 

-3655073397 

.4178090091 

.0263132347 

6 

.5421118538 

.5945791192 

.0196294939 

7 

.7601447375 

.7991787324 

.0130424499 

8 

1.0190503483 

1.0425922027 

.0080538340 

9 

1.3182897476 

1.3309553604 

.0046499131 

10 

1.6575617149 

1.6639218023 

.0025025914 

11 

2.0367717604 

2.0397792676 

| .0012503827 

12 

2.4559788735 

2.4573135565 

.0005779825 

13 

2.9153514543 

2.9159042425 

.0002465585 

14 

3.4151351711 

3.4153478013 

.0000968939 

15 

3.9556318607 

3.9557075370 

.0000350346 

16 

4.5371873692 

4.5372122223 

.0000116439 

17 

5.1601856915 

5.1601932080 

.0000035542 

18 

5.8250470008 

5.8250490909 

.0000009956 

19 

6.5322278278 

6.5322283615 

.0000002557 

20 

7.2822223566 

7.2822224815 

.0000000602 

21 

8.0755643170 

8.0755643438 

.0000000130 

22 

8.9128292551 

8.9128292604 

.0000000026 

23 

9.7946371142 

9.7946371151 

.0000000005 


n = 128 

n' = 32 


1 

.0055646485 

.0056474386 

.0000836954 

2 

.0286731117 

.0297432178 

.0011316178 

3 

.0686191075 

.0729961100 

.0047528127 

4 

.1238383560 

.1351282785 

.0115932142 

5 

.1941069901 

.2155251600 

.0183257845 

6 

.2813550805 

.3130880575 

.0204781205 

7 

. 3878324001 

.4262057173 

.0183965411 

8 

.5145492100 

.5533849095 

.0147608539 

9 

.6616147368 

.6949482493 

.0111662180 

10 

.8288591342 

.8538093471 

.0081292943 

11 

1.0160927402 

1.0331216673 

.0057288341 

12 

1.2231683248 

1.2341736061 

.0039093147 

13 

1.4499857472 

1.4568530792 

.0025797170 

14 

1.6964847002 

1.7006490926 

.0016434498 

15 

1.9626365278 

1.9650916323 

.0010092418 

16 

2.2484371391 

2.2498418722 

.0005966957 

17 

2.5539011953 

2.5546795382 

.0003393272 

18 

2.8790574811 

2.8794742689 

.0001854763 

19 

3.2239453612 

3.2241606823 

.0000973962 

20 

3.5886121959 

3.5887193728 

.0000491158 

21 

3.9731115618 

! 3.9731629082 

.0000237800 
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TABLE II.- Concluded 



n = 128 

n ' = 32 


22 

4.3775021010 

4.3775257587 

.0000110517 

23 

4.8018468266 

4.8018573035 

.0000049295 

24 

5.2462127401 

5.2462171976 

.0000021100 

25 

5.7106706464 

5.7106724677 

.0000008666 

26 

6.1952950885 

6.1952958030 

.0000003415 

27 

6.7001643509 

6.7001646199 

.0000001291 

28 

7.2253605012 

7.2253605984 

.0000000468 

29 

7.7709694533 

7.7709694870 

.0000000163 

30 

8.3370810429 

8.3370810542 

.0000000054 

31 

8.9237891120 

8.9237891156 

.0000000017 

32 

9.5311915985 

9.5311915996 

.0000000005 
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